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A Step by Step Computer Solution 
of Three Problems in Non-numerical Analysis 

This memo describes the step by step solution of three problems from 
different fields of applied mathematics. These problems are solved by 
typing a series of computer commands for the manipulation of symbolic 
mathematical expressions. These commands are best typed at the PDP-6 
console, so that the Type 30 display and the wider range of keyboard sym¬ 
bols can be used. The syntax of commands typed at the PDP-6 will be dev - 
scribed. These commands are translated into a string of symbols which are 
sent, to CTSS, where they are parsed into a LISP expression, which is then 
evaluated. 

The mathematical operators which are available in the system will be 
described and then the step by step solution of each of the problems will 
be given. 
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II. The Mathematical Operators 

The commands typed at the PDP-6 are similar to Algol statements. 

The commands are typed and executed one at a time. More complex operations 
involving the definition and alteration of commands and the introduction of 
more pneumonics and man machine interaction will be described later. The 
commands consist of infix operators, functions, and variables. Functions 
and variables can be subscripted and any subexpression can be quoted. A 
sample command which requires most of the notation is: 

//’E1«-'(X + Y) * ' DRV( ' :T, 1,DRV( 'U,2,El)) + !E2,20 + 1 (F[I,J] (X,Y))+2# 

In words: the name El is assigned to the expression which is the sum 
of three terms. The first term is the product of (X + Y) with the unevalu¬ 
ated first derivative with respect to lower case T of the second derivative 
with respect to U of the expression currently named El. The second term 
is the 20th subexpression of a displayed expression currently named E2; 
this subexpression has been indicated with the light pen. The final term 
is the square of a subscripted function of X and Y. The notation may seem 
somewhat complex, but as will be seen, a complex notation is required to 
express in a compact way the many small steps required to solve a particular 
problem. 
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The infix operators are: 

A+-B B is given name A. As such it is written 

on the disk. The value of +■ is B. 

!A,N The Nth subexpression of A is the value of !. Inten¬ 

sify the desired subexpression of A by pointing to 
its main connective with the light pen. Then type 
!A, and the computer will type N. If the expression 
has no main connective, point to one of its arguments 

and type ;!A instead of !A. Consider all minus signs 

... to be unary. 

A=B Equate A and B 

A+B A plus B 

A-B A minus B 

A*B A times B 

A/B A divided by B A/B*C is equivalent to A/(B*C) 

AfB A to the power B 

The functional, subscripting, and set notation is: 

A(C,D,E) A is a function with arguments C,D, and E. 

A[I,J] (C,D,E) A is a function with arguments C,D,and E. 

1, J 

A[I,J] A is a variable 

± , J 

(A,B,C) This is a set with three elements. By convention 


(A)= A 
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’ Either an expression or a variable name can be quoted 

A function name always stands for itself. Quoting 
a function name means that its arguments will be 
evaluated but that the function will not be evaluated 
For example let F(X,Y) = X-Y be a function and let 
X and Y be names for A; then 
F(X,Y) evaluates to 0 
'F(X,Y) evaluates to F(A,A) 

F(’X,Y) evaluates to X-A 
'(F(X,Y)) evaluates to F(X,Y). 

Quoting a function or variable name does not quote 
its subscripts. Numbers are taken as quoted auto¬ 
matically . 

: Causes the letters which follow it to be lower case 

for purposes of display. 

As in CTSS, there are two editing characters: 

? Deletes all the characters of a command back to the 

initial #. 

Deletes only the immediately preceding character. 

"//” must be the first and last character of every command. "; n 
causes the current intensified subexpression to be raised one level. For 

B R 

example, if the A in A + C is intensified, then when ; is typed A will 


be intensified. 
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Other available 
ALLSUMEXPAND(EXP) 

BRINGOVER(EXP, X) 

COLLECT(EXP,SET) 

DEPENDENCE(EXP) 

DELSUBST(EXP,OLDDEL, 

DRV[X ,N , ...» X ,N 
11 n i 

DRVDO(EXP,X) 

DRVFACTOR(EXP, X, N) 

DRVZERO(EXP,X) 

EVALUATE(EXP,SET) 


operators are: 

Applies SUMEXPAND to every summation in expression 
EXP. 

Subexpression X, which has been indicated with the 
light pen is brought to the other side of equation 
EXP. 

Top level terms in EXP are collected on powers of 
the expressions in set SET. 

Returns a set of the variable and function names in EXP. 


NEWDEL) dx_ ^ dx_ 

d OLDDEL ^ d NEWDEL 

subexpression in EXP 


for each such 


, Y) Differentiate Y times with respect 

to X., for each i. 

1 

All indicated derivatives with respect to X in EXP 
are carried out as far as possible. 


d N+M f d M d N 

—££ ( —) for each such subexpression in 
dx dx dx 

EXP. 

All derivatives with repect to X in EXP are set 
equal to zero. 

SET is a set of equations; whenever the left side 
of one of these equations can be matched to a sub¬ 
expression in EXP, the right hand side is substituted. 
The left sides must be variables or functions. 

A match occurs whenever a binding of the function 
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r 

\. 

EXCHANGE(EXP) 

EXPAND(EXP) 

FACTOROUT(EXP, FACTOR, 

GROUP(SET) 

S' 

V 

LEFT(EXP) 

LIMIT(EXP,X,N) 
MULTIPLYTHROUGH(EXP,X 
NEWNAMEO 
NORMPOLY(EXP,X) 

REPLACE(E,X,Y) 


variables and subscripts can be made. 

If the top level connective of EXP is binary, its 

arguments are exchanged, right to left. 

Multiplies out all expressions of the form a*(b+c) 

in EXP. In addition, 

d / ., N da , db 
— (a+b) -> — + — 

dx dx dx 

Y) The factor FACTOR is factored from each term of 
EXP. The third argument Y is optional. If Y is 
present, the factor FACTOR is renamed Y. 

The set SET of terms which have been indicated by the 
light pen in EXP are grouped within the associated 
sum or product. The value of GROUP is the grouped 
set of terms. 

Returns the left argument of the main binary connec¬ 
tive of EXP. 

Determines the limiting value of EXP as X approaches N. 

) Multiplies each top level term of EXP by X. 

Creates a name of the form Fn, where n is an integer. 

Every sum in EXP is treated as a polynomial in X 
and a power of X is factored out so that the lowest 
power of X in the polynomial will be zero. 

Expression X replaces Y in the expression named E. 


Y is a term indicated with the light pen or a group of 
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tenns indicated with GROUP. If the light pen has 
been used to construct X, the resulting expression 
position is named HOLE. HOLE can then be used for 
the third argument. If X is equal to NIL, then the 
third argument is omitted from the expression named E. 
RIGHT(EXP) Returns the right argument of the main binary connec- 

tive o f EXP. 

SIMPLIFY(EXP) Simplifies expression EXP. 

SOLVE(EXP,X) Solves equation EXP for variable X as far as possible. 

SPLIT(EXP) Subparts of EXP are named and replaced by their names 

in EXP, so that EXP will contain less than 100 sub¬ 
expressions . 

SUBSTITUTE(EXP,X,Y) Substitute X for each occurrence of Y in EXP. 

SUMEACH (EXP) Z(a+b) la + Zb 

SUMEXPAND(EXP) Expands the finite summation EXP. 

TERM(EXP,N) Returns the Nth argument of the top level connective 

of EXP, or NIL if there is no Nth argument. 

TRUNCATE(EXP,VAR,N) Expands EXP up to power N in variable VAR. 

SUM(I,N1,N2,Y) Sum expression Y for values of I from N1 to N2. 

ITG(X,L1,L2 ,Y) Integrate Y with respect to X between limits LI and L2. 

Expressions which are assigned names are kept on the disk. The expres¬ 
sion most recently computed always has the name LAST. When A«-B is executed, 
if A is not ,! LAST'' and is already the name of an expression, then this old 
value of A is given the name OLD. Thus, if A+A+2 is executed and then is 
found to be incorrect, the old value of A can be retrieved. 



Operators used for input-output and disk storage are: 


EDISPLAY(E) 

EPRINT(E) 


EDELETE(E) 


Displays the expression named E on the PDP-6 scope. 

Prints out the internal form of the expression 
named E with PLS, PRD, EQN, and PWR in infix form; 
the other operators in prefix form. 

Deletes expression named E from the disk. 


This completes the description of the PDP-6 commands. 




The Poincare-Lighthill Procedure 

2 3 

Applied to x + 10 x = lx 


The Poincare-Lighthill procedure is typical of a number of procedures 
used to find the first few terms in the asymptotic expansion of the function 
which is the solution to a mildly non-linear differential equation. The 
equation chosen here is that for a harmonic oscillator with a small forcing 
function. These solution procedures involve assuming a series expansion 
in powers of the small parameter for one or more of the parameters and 
variables in the equation, substituting these series into the differential 
equation, and thus obtaining a series of relations between the coefficients 
of like powers of V. Each of these equations is then treated in turn by 
whatever methods seem appropriate. Thus it is in general necessary to see 
these equations before the next steps in the solution process can be deter¬ 
mined. 


When a typed command has been completed, the machine makes a response 
of acknowledgment. This standard response will be omitted in the dialogue 
to follow; only the typed commands and the displayed equations will be shown. 
A running discussion of the dialogue is included and the displayed equations 
are shown on the pages following their use. These equations were plotted 
by the CALCOM? plotter. A photograph of these same equations is shown at 


the end 


section. 


reader should be 


that 


equation 


syntax used, more than 


line 


expres; 


occur over a divide 


witnn 


T n i s 


i 1 lust rat ted 


ecua.tr on 


sect ion. 
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Enter the differential equation. 

//’El«-' (DRV(:T,2,X(:T)> + OMEGA +2*X( :T) = EP*(X( :T)) +3)// 

//EDI SPLAY (’El)// 

A new independent variable x is introduced in order to stretch the 
time. Type in expressions for series expansions for X and t in terms of 
functions of t. That is, a solution for X(t) rather than for X(t) will 
be found. Since t depends on x, equation El can be used to find an equa¬ 
tion in derivatives of X(x) . As a final step the inverse relationship 
x(t) will be found, so that X(x) will give X(t). 

#’E2-<-' (X(TAU) = SUM(I,0,INF,EPtI*X[Ij (TAU)))# 

//EDISPLAY ( ' E2) # 

//*E3«-' ( :T(TAU) = TAU + SUM(J ,1, INF,EP+J* :T [ J] (TAU) ) )// 

//EDISPLAY (’E3)// 

In order to substitute dx for dt it is necessary to apply the trans¬ 
formation 

d 2 x d_ ^ dx ^ 

dt 2 dt dt 

to equation El. 


//’E4^DRVFACTOR(El, ' :T,1)// 
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Display E4 for comparison with the substituted result below. 
//EDISPLAY (£E4) // 


Now substitute dx for dt, and X(x) for X(t) 


#'E5^SUBSTITUTE(DELSUBST(E4 , 1 (DEL(:T) ) , ’(DEL(TAU))/DRV('TAU,1,RIGHT(E3))), 

’(X(TAU)) ,' (X( :T) ) ) // 


//EDISPLAY('E5)# 


Now substitute the series for X(t) and perform the indicated differentiation 
with repect to x. 


//'E6^DRVD0( SUBSTITUTE (E5, RIGHT(E2),'(X(TAU))),’TAU)# 
//EDISPLAY ( f E6)// 


Now expand both sides to first order in e. 

// f E7^TRUNCATE (E6, ' EP, 1) # 

The zero order terms form the harmonic oscillator equation; the solution 
can be written down by inspection as Acos cot. Use the light pen to form 

an equation of the first order terms. 

#'E8^!E7,6=;!E7,88# 

//EDISPLAY (’E8)# 

Bring the terms in X i (x) to the left side of the first order equation. 

Substitute for X (x) and carry out the indicated differentiation. 

0 



(})x- s =(l)x- (^)x 

E 



vr 






# f E9+-SIMPLIFY(DRVD0(SUBSTITUTE(SOLVE(E8, ’ (X[l] (TAU))) ,' (A*COS(OMEGA*TAU)) , 


* (X[ 0] (TAU) ) ) , 'TAU) ) # 


//ED I SPLAY (' E9) // 

It is now necessary to substitute an identity for cos 3 a)x and to collect 
terms on sinwx and cosmx. 

//'El(KCOLLECT(EXEAND(SUBSTITUTE(E9, f ((COS(3*OMEGA*TAU)+3*COS(OMEGA*TAU))/4) , 

’ ( (COS(OMEGA*TAU) ) f3) ) ) , ' ((SIN(OMEGA*TAU) ,COS(OMEGA*TAU) ) ) ) // 

//EDISPLAY (' E10) // 

Theoretical considerations require that 

the coefficients of coscox and sinmx must be zero if there is to be a 

periodic solution for X (x). From the coefficient of sinwx it is apparent 

that t'(x) must be some constant C. This constant is determined from the 
1 

coefficient of coswx. 

//'Ell-<-SIMPLIFY(SOLVE(SUBSTITUTE(!E10,44, 1 C, ' (DRV(TAU, 1, :T[1] (TAU))))=0, ’C))// 

//EDISPLAY('JSU)# 

3A 2 

So X = A coswx and to first order t = x + e- x. 

0 8 co 2 

Thus to first order x « t (1 - - ) and X A cos (1 - - )t 

8w 2 0 8o) 2 
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One effect of the nonlinear term is thus seen to be an alteration in 
the frequency of the zero order term. 

In conclusion, note the large number of small steps necessary to 
solve this problem. These are the result of doing almost the entire solu¬ 
tion in the machine. In the case of a small problem such as this, some 
of these steps could be done by hand, the object here is to illustrate the 
steps which would be required for a larger problem. This problem also illus 
trates how rather lengthy intermediate calculations can lead to some rather 
concise results. The perturbation of the frequency in the zero order func¬ 
tion is found to have a simple expression. 










f 
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III. Plasma Accelerator Electrode Boundary Layers 

The second problem is a duplication of the work in the first three 
sections of chapter three of a 1963 Masters Thesis by J.S. Draper for 
the M.I.T. Department of Aeronautics and Astronautics. This thesis inves¬ 
tigates the laminar compressible boundary layer on the electrode walls 
of a direct-current crossed field plasma accelerator under very special 
physical conditions. Many of the assumptions used under these conditions 
are set forth in the paper "Electrode Boundary Layers in Direct-Current 
Plasma Accelerators" by Jack L. Kerrebrock in the August 1961 issue of 
the Journal of Aerospace Sciences . Kerrebrock's paper investigates a 
solution involving less mathematical manipulation than that undertaken by 
Masters Thesis student Draper. 

In summary, the entire solution procedure is as follows: 

1. Write down 5 non-linear partial differential equations: 

Momentum 

State 

Continuity 

Energy 

Electron mobility as a function of temperature. 

These equations relate U stream velocity 

V lateral velocity 
t temperature 
p density 

p electron mobility 


P pressure 
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in terms of the independent variables x and y. The constants are: 

j current 
B magnetic field 
Cp specific heat 
K compressibility 
G conductivity 
R gas constant. 

The absence of variation in the y direction in the free stream is used 
to find the momentum and state equations there. These two reduced 
equations are solved for — which is eliminated from the 5 main equa¬ 
tions, since P is not a function of y. 

The relation H = C t + U 2 /2 is used to substitute derivatives of H for 
those of t in the energy equation which then becomes an enthalpy equa¬ 
tion. t is thus eliminated from this equation. Simplification of the 
resulting expression requires introduction of the momentum relation. This 

step is performed because the enthalpy equation has a term proportional 
± C y 

to [1 - p ] where P^. is and can be approximated as 1, thus elimina- 

r ■ 

ting this term. 

Next, a change of independent and dependent variables is made. The change 
of independent variables is such that it approximates a similarity trans¬ 
formation for low Mach number. These transformations change x and y to 


x and n. In addition 


t' H 

is defined as f 1 , — as 0, and — as e 

t H 6 


The momentum and enthalpy equations are transformed, using the continuity 


and state equations as side conditions, x is then changed to 


There 


result two non-linear differential equations in f and g and their deri¬ 
vatives with respect to n and M 
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\ 

5. f and g are then approximated as f ? (n,M ) = b(M )r\ + c(M')n 2 , 

00 00 00 

g(n,M ) = e(M)n + f(M )n 2 . These approximations are substituted 

00 00 00 

into the two non-linear differential equations. The equations are 
then integrated with respect to n between the wall and the edge of the 
velocity boundary layer 6u and the edge of the enthropy boundary layer 
<5e respectively. bCM^) and eCM^) are eliminated from the result by 
the relations f’CSujM^) = 1 and gC&ejM^) = 1. There result two ordin¬ 
ary linear differential equations for the derivatives of 6,u, c, &e, 
and f with respect to M^. 

6. Two more linear differential equations for 5u, c, 6e, and f are generated 
by choosing the coefficients of the approximations in step 5 so as to 
satisfy the momentum and enthalpy equations produced in step 4 exactly 

at the extremal points f" = 0 and g’ =0. 

7. The four resulting linear differential equations are solved for the 

96u 3c 36e , 3f 

derivatives—by Gaussian reduction. These 

00 00 00 00 

four expressions are then numerically integrated with a Runge-Kutta 
method. 

This problem has several interesting features. It is a demonstration 
of the notation use by workers in this area. The algebraic expressions 
are of a size difficult to manipulate by hand, but within the capabilities 

of current machines. The final symbolic result is large; it is difficult to 
write the corresponding numerical integration program correctly when this 
result must be input by hand, but here it is developed in the machine and 
could then be transformed into the required numerical program. Note that 
the symbolic steps are needed in order to cast the problem in terms of the 
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independent and dependent variables of interest. The problem is charac¬ 
terized by the application of simplifying side conditions and physical 
assumptions. As such, it involves a number of manipulations for the pur¬ 
pose of expression condensation. This will be apparent from the following 
step by step reproduction of the first three sections of Chapter III. 

These steps bring the described solution through the application of the 
similarity transformation to the momentum equation. 

Input the momentum equation: 

#’Dl^ f (RH0*(U*DRV(X,1,U) + V*DRV(Y,1,U)) = DRV(Y,1,MU* 

DRV(Y,1,U)) - DRV(X,l,P) + :J*B)# 

Input the energy equation: 

# , D2«-'(RH0*C[P]*(U*DRV(X,1,:T) + V*DRV(Y,1, :T)) = 

DRV (Y,1,K*DRV(Y,1,:T)) + MU*(DRV(Y,1,U) ) +2 
+ U*DRV (X, 1, P ) + :Jt2/SIGMA# 

The boundary layer solutions must match the free stream solution. 

The free stream values are indicated by the subscript °°. At the free 
Stream, there is no variation in the boundary layer with respect to y. 

To save rewriting, define sets containing the variables to be subscripted. 
#'D3^'(RHO,SIGMA,U,:T,H)# 

#'D4-*-' (RHO [ INF ] , SIGMA [INF] ,U[INF] , :T[INF] ,H[INF])// 

Then in the free stream D1 and D2 become: 

#'D5^SIMPLIFY(SUBSTITUTE(DRVZERO(Dl,’Y) ,D4,D3))// 

# ’ D6-*-SIMPLIFY ( SUBSTITUTE (DRVZERO (D2 , ’Y) ,D4,D3))// 

#EDISPLAY(’D5)r 


/REDISPLAY ( 'D6)/- ; 







D5 and D6 will now be used to express some of the expressions in D1 and 

D2 in terms of the free stream quantities. In order to substitute 

for some of the terms in a sum, the terms must be grouped using the light pen. 

This is somewhat inconvenient. 

//EDI SPLAY (Dl) # 

Now the last two terms in Dl are replaced by the left side of D5. 
//'D7+-REPLACE ('Dl,LEFT(D5) , GROUP( (!Dl,33,!Dl,38) ) )// 

#EDISPLAY('D7M 

#'D8-*-LEFT(D2) = EXPAND (SUBSTITUTE ( 

RIGHT (D2) , RIGHT ( SOLVE (D6,' (DRV(X,1,P) ) ) ) , ’ (DRV(X,1,P))))// 

//EDISPLAY ('D8) // 

The last two terms in D8 are now put in a factored form. 

#'D8+-REPLACE (' D8, FACTOROUT (GROUP ((! D8,62,! D8,78) ) , 1 ( i J+2/SIGMA[INF] ) ) 

' HOLE)# 

#EDISPLAY('D8)# 

In order to effect a cancellation, equation D8 will now be transformed 

by replacing the temperature, t, with the enthalpy, H, using the definition 

H = Ct + U 2 /2. 

•P 

First H = C p t + U 2 /2 is solved for temperature t, 

#'D9«-S0LVE(' (H(X,Y) = C[P]*:T(X,Y) + U(X,Y)t2/2) , ' ( :T(X,Y)))# 

Now the substitution is made for t and t . 

oo 

# ' DIO^-EXPAND (SUBSTITUTE(D8 , (DRV ( ? X, 1, RIGHT (D9) ) , SUBSTITUTE (DRV ( ' X, 1, RIGHT(D9 ) ) , 
D4,D3),DRV('Y,l,RIGHT(D9))),’( DRV(X,1,:T), DRV(X,1,:T[INF] ) , 

DRV(Y,1,:T))))# 



W/j„. 


■ *?«#**•> 


•• * ■ *- 4 Vw»*6- - W j •■ **» >Si.'-*' ft,! T.'> ' •-*« 


•-i.-"^ i • ■ ‘ft**** <-**■• -»• • ^ «»* < i frw fttur nrw i l lii ti i 




■ '. r ir^-.4 ■'■■’• ‘"**3 


**■ 


3 > 
D D 


1 —3) b b b 
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^ b’ 


=3 3 b 


3 > 
U D 


ID > 

U U 


3 > 
U D 


<-> > 
V U 


*-> > 

u u 


3 X 

u u 


^ > 
u z 


4 -> > 

TJ U 


+-> X 
U D 


4-> X 
U D 


/“V> 
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Now, there are a number of tedious grouping steps. 

//EDISPLAY( ’DIO) // 

First, factor RHO out of two of the terms on the left side of the equation 
// 'D11+-REPLACE ( ’DIO, FACTOROUT (GROUP ( ( ! DIO ,48J DIO, 21))~ T RE0) , ’HOLE) // 
//EDISPLAY (’Dll)# 


Now factor RHO*U from the other two terms and bring them to the right side 
#' D11+BRING0VER( ’Dll,FACTOROUT (GROUP ( ( ! Dll, 5,! Dll, 45) ) , ’ (RHO*U (X,Y) ) ) ) # 


A machine matching operation would be better for the next factoring 
step which must be done twice when light-pen pointing is used. 

//EDISPLAY (’Dll)# 


The quantity K/C p P is factored out of two of the terms and set equal 
to 1/P r . 

//'D11+-REPLACE( 'Dll,FACTOROUT(!Dll,34, ’ (K/ (C [P]*MU)) , ' (1/P[R])) , ’HOLE)# 
//EDISPLAY (-’Dll)// 

//’Dll-HREPLACE( 'Dll,FACTOROUT(!Dll, (K/(C[P]*MU)) , ' (1/P[R])) , ’HOLE)// 

C^p/K is the Prandtl number P r * It is close to 1, and setting it 
to 1 will effect a simplification if the identity 


/^ U \ 2 

^ = 


dY 


[ U .ucx.y^LII] - 


U 


d r dU, 

d7 [v d? ] 


is also substituted. This simplification means that the heat conduction 
away from a point is just equal to the viscous dissipation at that point. 
#'D11^SIMPLIFY(SUBSTITUTE(Dll,’(1,DRV(Y ,1,MU*U(X,Y)*DRV(Y,1, 

U(X,Y) ) ) -U(X,Y)*DRV(Y,1,MU*DRV(Y,1,U) )) , ' (P[R] ,MU* (DRV(Y ,1 ,U) ) + 2) ) ) // 
//EDISPLAY( ’Dll)# 

Next, the substitution in Dll of the right side of equation D7 for 


its left side effects a nice cancellation. 









u(x, v) *nu 










//’Dll-f-SIMPLIFY(REPLACE( ’Dll, (-RIGHT(D7) ) ,GROUP( (!Dll, 117,!Dll, 122) ) ) ) # 

//'D11*TjEFT(D ll) = SIHPLIEY(EXPAND(SUBSTITUTE(RIGHT(Dll) , ' (U,U[INF]) , 1 (U, (X,Y) , 

U[INF] (X,Y) ))))// 

//EDISPLAY(’Dll)// 

The momentum equation D7 and the enthalpy equation Dll are now 

ready for the similarity transformation. 

The transfoomationsodf independent variables are: 

//'D12«-’(XB(X) = ITG(X,0,X,P*U[INF]/ (P[0]*U[0] ) ) )// 

#’D13^’(ETA(X,Y) = (U[INF]/U[0])*((U[0]/(2*NU[0]*XB(X)) + (l/2))* 

ITG (Y ,0, Y ,RHO/RHO [0 ]) )// 

Next, to compute the required differentials. 

First 1 = dn . JL ] dX . _ 1 

dX dX dn + dX dX 

//'D14*rDRV ('X, 1,RIGHT (D13) ) *' (DEL (ETA) )+DRV(’X; 1,RIGHT (D12) ) * ’DEL (XB?) ) // 
//EDISPLAY ( * D14 ) // 

D14 now contains the expression for ryas a factor. 
//*D14+REPLACE('D14,FACTOROUT(!D14,4,RIGHT(D13) , ’ (ETA(X,Y))) , ’HOLE)// 

//EDISPLAY(’D14)// 

n iff now substituted for its definition and (P*U oo )/(P 0 * U 0 ) is 
factored out. 

//'D14+-SIMPLIFY(FACTOROUT (SUBSTITUTE (D14,DRV ( 'X, 1,RIGHT (D12) ) , 

' (DRV(X,l ,XB (X) ) ) ) , 'P*U[INF] / (P [0] *U[0] ) ) ) ) // 

//EDISPLAY (’D14)// 

The differential for 1/dY is 
// ’ D15+DRV( ’ Y, 1, RIGHT (D13) ) * ' (DEL (ETA) ) // 












#EDISPLAY(’D15)# 


Expressions D14 and D15 are the required differentials for ~r and 

QA dY 

respectively. 

Now to transform the momentum equation D7. First, substitute for 
the differentials of the independent variables and the new normalized 
dependent variable f, defined by U/U^ = f'CXin). 
#'D16^SIMPLIFY(SUBSTITUTE(DELSUBST(DELSUBST(D7,*(DEL(X)),D14), 

'(DEL(Y)),D15),*(U[INF](XB)*DRV(ETA,1, 

:F(XB;ETA)) ,U[INF] (XB) ,XB,ETA) , ' (U,U[INF] ,XB(X) ,ETA(X,Y))))# 

Now assume y= y t/t . 

0 0 

#'D17^SIMPLIFY(SUBSTITUTE(’D16, ? (HU[0]*:T/:T[0],’MU))# 


Use the equation of state to recognize that 


k f ' 0 


since P is taken independent of Y, 


beemade. 


A substitution for P should therefore 


#EDISPLAY('D17)# 

The machine responds with ’’EXPRESSION TOO LARGE" 

Therefore, the left side of equation D17 will be treated first, the 
substitution for P will be deferred. 

#' Dl8-«-LEFT (D17 ) # 

An expression for V will now be developed and substituted into D18. 

It was shown in Chapter II of the thesis that the continuity equation yields 



i) ■ eta(x, v)- del(eta) 

- Mxj -™ l(xb) 




Enter this expression. 

// t D19^ , (-(RHO[Q3/RHO)*DRV(X,l,((2*U[0]*NU[O]*XBf(l/2))*:F(XB.ETA)))# 

Now substitute the new independent variables: 

# , D20+-SUBSTITUTE(DELSUBST(D19, 1 (DEL(X)),D14), 1 (U[INF](XB),XB,ETA), 

*(U[INF],XB(X),ETA(X, Y ) ) ) # 

Now, substituting this expression for Vninto the partially transformed 
left side of the momentum equation, D18, the entire expression is differen¬ 
tiated as far as possible, expanded and simplified. 

#'D21^EXPAND(DRVD0(DRVD0(SUBSTITUTE(D18,D20, 'V),’ETA),’XB))# 

Now some shorthand will be introduced to make D21 easier to read. 

#'D22-*- 1 ( : F,U, FI,F2, F01,Fll,U[ INF] ) 

# , D23^' (:F(XB,ETA) ,U(XB) ,DRV(ETA,1, :F) ,DRV(ETA,2, :F) ,DRV(XB,1, :F) ,DRV 
(ETA,1,XB,1,:F),U[INF](XB))# 

# , D21^-SUBSTITUTE(D21 ,D22 ,D23) # 

#EDISPLAY('D21)# 

It is now convenient to factor a large coefficient from each term 
and set it equal to 1 since it will also be factored from the other side 
and set to 1 there. 

#'D24*-SIMPLIFY(FACTOROUT(D21, f (RHO*U[INF]+3*P/(2*P[0]*U[0]*XB)),1))# 

#EDISPLAY(*D24)# 

Returning to the right side of D17. 

#'D25+-RIGHT(D17)# 

#EDISPLAY('D25)# 

The deferred substitution for P is now made. 
#'D27^REPLACE(’D25,FACTOROUT(!D25,44, f ( (RHO* : T)/(RHO[0]*:T[0])),'(P/P[0])),’HOLE)# 
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The expression is now differentiated as far as possible,as was the 
left side. 

# ’ D28-HEXPAND (DRVDO (DRVDO (D27 , ' ETA) ,' XB) ) # 

The simplifying notation is substituted. 

// 'D28+-SUBSTITUTE(D28 ,D22 ,D23) // 

REDISPLAY ( , D28)// 

Finally the large factor is removed from each term as it=was from 
the right side. 

# ’ D2 9+SIMPLIFY (FACTOROUT (D2 8,’ ( RHO*U [ INF ] + 3*P / ( 2*P [ 0 3 *U [ 0 ] *XB ) ) , 1) ) // 

The two sides of the transformed momentum equation are now recombined. 

# ’D30-*-S IMPLIFY(D24-D29) = 0# 

//EDISPLAY (’D30) # 

Using some more relations developed earlier in the thesis, a final 
substitution will be made, from the equation of state at constant pressure. 

P X. i. P 

oo l U 00 

-— = r— but —— =0 so — = 0 , and from Chapter II, 


0 = (0 - 0 )g + © - (0 - l)f 2 

s w ° w s 


The final transformed equation is: 

# ’D31+SIMPLIFY (SUBSTITUTE (REPLACE ( ’D30, FACTOROUT (! D30,38 , ’ (RHO [ INF] /RHO) , 

’ ( (THETA[S] - THETA[W] )*G + THETA[W] - (THETA[S] - l)*Fl + 2) ) ,’HOLE) , 
’ (1 + (GAMMA - l)*M[INF]t2/2) , ’ (THETA[S] ) ))# 


#EDISPLAY(’D31)// 






In conclusion, since physical approximations are involved in the 
expression condensation, close man-machine interaction is required. The 
greatest draw backs to sufficient interaction are the input notation and 
the lack of sufficient facilities for abbreviation on output. 



IV. A Multiple Channel Queueing Problem 


The third problem is taken from Interim Technical Report No. 13 of 
the M.I.T. Center for Operations Research titled "A Class of Queueing 
Problems". This work is a 1955 Doctor’s Thesis by H.N. Garber. The 
queueing situation shown by the example in Figure 1 is treated. 


arrivals 



4 

3 

2 


4 

• 

3 

2 

i 

4 

3 


i 


Channel 1 


Channel 2 


Channel 3 


Arrival times at the queueing complex are exponentially distributed with 
mean A. Arrivals enter any channel which becomes vacant and progress 
through several phases of service with exponentially distributed service 
times with mean ky.' There are M channels and k 'phases of service in the 
general case. It is desired to find the probability distributionnof the 
number of units in the system. The solution can be carried a number of 
steps for the general case. A set of equations relating the probabilities 
of the states are written, where a state of the system is taken as a cer¬ 
tain number N in the system and a description of which phases of which 
channels are occupied. Each equation is then multiplied by the appropri¬ 
ate variables so that the summation of the equations will yield an equa¬ 
tion with the generating function for the state probabilities as its left 







side. Unfortunately, the generating function also appears in a summation 
on the right side. The independent variables of the generating function 
are constrained to make each term of this summation equal to zero. This 
constraint is expressed in a change of independent variables. It is then 
shown that a summation over values of the new independent variables will 
yield the old generating function. 

This generating function still involves the state probabilities for 
a partially full system as eonstantsttotbe determined. There are a number 
of relations which can be used to determine these initial probabilities. 

As the third problem this evaluation will be carried out for the case of- 
two channels with two phases of service. 

The source of expression complexity is different in this problem. In 
the first problem complexity arose because a small parameter allowed an 
unsolvable problem to be split into a spectrum of solvable ones. In the 
second problem complexity was the result of the number of important pheno¬ 
mena in the physical process being described. Here, complexity is the 
result of the number of states in the process being described, all of 
which are very similar. In this problem there would be the best chance 
for reduction in complexity through proper notation. 

This problem has been included so that the reader can better refine 
his intuitive notion of the types of mathematical operations and notation 
needed to solve problems in different areas of applied mathematics. Sum¬ 
mation expansion, function evaluation, limits, and some more grouping 


operations are introduced. 




The function in the transformed variables is: 

(K(Z,Q[1],Q[2]) = ((Zf(2*:K + 1))*SUM(R,1,:K,(ZfR)* 

A[R] * (:E+(2*PI*I*Q[1]*R/ ;K)+:E+(2*PI*I*Q[2]*R/ ;K))))/ 

(2*(Z*W(Z+:K) - ALPHA(Q[1],Q[2]))))# 

Where W, a, and A are defined by: 

r 

#'Q2^ f (W(Z+:K) = 1 + THETA - THETA*Zf :K)// 

# , Q3^ , (ALPHA(Q[1],Q[2J) = (:E+(-2*PI*I*Q[l]/:K)+ 

:E+(-2*PI*I*Q[2]/:K))/2)# 

(A[R](Z) = P(1,R+1,0) - (2*W(Zt:K) - 1)*P(1,R,0))# 

Redisplay (' qi ) # 

It is useful to have the denominator of Q1 contain only powers of Z . 

k-1 k-1 k _ k 

This is done by using the identity (X-Y)(X ...Y ) = (X - Y ). At 
present the system contains no operators for achieving this goal, so it 
must be done by brute force. 

#'Q5-*-LEFT(Ql) = ' (SUM(J,0, :K-1, (Z*W(Z+ :K) ) +J*ALPHA(Q[1] ,Q[2] f ( :K-1-J) ) )* 

SUBSTITUTE ( RIGHT (Ql) , ’ ( (Z +: K) * (W ( Z t: K)+ ; : K) - ALPHA (Q[l],Q[2])t:K) 
' (Z*W(Zt :K) - ALPHA(Q[1],Q[2])))// 

#EDISPLAY('Q5)* 

Now substitute 2 for the number of phases of service, Then 

expand the summations. 

#'Q6^SIMPLIFY(ALLSUMEXPAND(SIMPLIFY(SUBSTITUTE(Q5,2, ’ :K))))// 
y/EDISPLAY ( ' Q6) // 



-0, *S 2*<n*I*0 



6 - 0-2 +1 







The substitution - lei — 2 is made in all the initial equations as well 

I 

since theywwlllbbe used several times. 

#' Ql+SIMPLIEY(SUBSTITUTE(Ql,2, t :K))» j 

# , Q2^ SIMPLIFY(SUBSTITUTE(Q2,2:K))# 

#’Q^SIMPLIFY(SUBSTITUTE(Q3,2,' :K))# 

# ’ Q4fSIMPLIFY ( SUBST ITUT E(Q4 ,2, 1 :K) )# " 

Now write down an expression for the generating function as a summa¬ 
tion of Q6 over the transformed independent variables, the q’s. 
#3Q7^(H(Z+2,V[1],V[2]) = (1/4)*SUM(Q[1],0,1, 

SUM(Q[2],0,1, ((V[l]/Z)*(-1)+Q[l] + V[l]+2/(Zf2))* 
((V[2]/Zl*(-l)tQ[2] + V[2]t2/(Zt2))*K(Z,Q[l] # 9[2]))»# 

Redisplay(’Q7)# 

Now expanding the summations, substituting the appropriate values 
of Q6, Q3, and and simplifying. 

#$Q8*-SIMPLIFY(EVALUATE(ALLSUMEXEAND(Q7) , (Q6,Q3) ) ) // 

REDISPLAY('Q8># 

To evaluate p(1,1,0) and P(1,2,0) certain known conditions are next 
imposed. Q8 is known to have a zero of order four in Z for all values 
of Vi and V2, so one would like to collect terms on Z 4 . The simplest way 
to explore this would be to expand Q8 and collect terms on Z 4 . Unfortunately, 
this leads to roughly a sixteen~fold growth in expression size and to memory 
overflow. 

Inspection of Q8 shows that it might be rearranged while in factored 


form. As the first step, subexpressions which are polynomials in Z, Vj, or 
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V 2 , have factors of Z, Vj_, or V 2 removed so that their lowest order term 
is of zero order in these. 

// 'Q81-«-SIMPLIFY(NORMPOLY(NORMPOLY(NORMPOLY(Q8 , ’ Z) , ’ (V[l] ) ) , 1 (V[2] ))) // 

//EDI SPLAY (’Q81)// 

Next the center two terms are combined. 

//' Q82-HR.EPLACE ( ’ Q81, FACTOROUT (GROUP ( (! Q81,85,! Q81,108) ) , 

’ (A[2]/W(Z+2))) , ’HOLE)// 

//ED I SPLAY ( ' Q82) # 

The resulting term is then expanded. 

//'Q83-«-REPLACE('Q82, EXPAND(! Q82,90) ’HOLE)// 

//ED I SPLAY ( ’Q83)// 

Now the other two terms are combined. 

//' Q84+REPLACE ('Q83, FACTOROUT (GROUP ( (! Q83, 111,! Q83,30) ) , 1/Q83! 70) , ' HOLE) // 
//EDISPLAY(*Q84)// EXPRESSION TOO LARGE 

Q84 will not display, so it is reduced in size by naming a subpart. 

//' Q85+-SPLIT (Q84 ) // 

//EDI SPLAY ('Q8 5)// 

Next the numerator of the larger term in 085 is arranged on powers of Z. 

// ’ Q86-<-REPLACE( ’Q85, COLLECT (EXPAND (SUBSTITUTE (! Q85,31 ,F2, 'F2) ) , »Z) , ’HOLE)# 
//EDISPLAY( ’Q86) // (EXPRESSION TOO LARGE) 

//’087^SPLIT(Q86)// 

//EDI SPLAY (’Q87)# 

Forming a term from the renamed parts of Q87 one has for the coeffi¬ 
cient of Z 2 in 087. 

//'088^SIMPLIFY(F2 + F3 + F4 + ' (2*Ai) + F5 + F6)# 
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#EDISPLAY(’Q88)# 

Looking at equation Q87, one can see that in a double power series 
expansion in Vi and V 2 , only the coefficient of V 2 V 2 does not have a 
zero of order 4 at Z = 0. Setting this coefficient equal to zero one 
obtains by inspection of Q87: 

A 2 Aj + Z 2 • W(Z 2 >A 2 

2W(Z 2 ) + Z'(Z 2 • W z (Z 2 ) -1) = ° 

Enter the equation into the machine. 

#'Q9+- T (A[2]/(2*W(Z+2)) + (A[l] + Zt2*W(Z+2)*A[2 J)/ 

(2*(Z+2*(W(Zf2))f2 -1)) = 0)# 

Evaluating W, Ai and A 2 one obtains at Z = 0 
# ' Q10^-SIMPLIFY (SUBSTITUTE (EVALUATE (Q9,(Q4,Q2)),0,'Z))# 

#EDISPLAY ( ’ Q10) // 

Recognizing that by definition P.(l,3,0) = P(1,1,0) 
this equation can be solved for P(1,2,0). 

#'Q11*-SIMPLIFY(SOLVE(SUBSTITUTE(Q10, ? (P (1,1,0) ) , » (P (1,3,0) ) ) , 1 (P (1,2,0)))) # 
#EDISPLAY( 1 Qll)# 

011 can be simplified somewhat. 

//' Qll-^-LEFT (Qll) = (-SIMPLIFY (FACT0R0UT(! Qll, 11, ! Qll, 29)/ 

FACT0R0UT(!Qll,43,!Qll,49)))# 


//EDISPLAY( 'Qll)// 




(i+e) (i+e) 
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2 

Let u « z , the generating function for the unconditional state 

00 

probabilities G(u) = u n p(n) can now be written. There are special 

n=0 

case terms for n = 1, and n = 2, The other terms are found from H(u,l;l). 
The most compact formula for H is Q86. 

//'Q12+-'(G(U)) = '(P(0,0,0)) + 2*’U*(' (P(1,1,0)) + RIGHT(Qll)) 

+ SIMPLIFY(SUBSTITUTE(EVALUATE(RIGHT(Q86), (Q2, Q4)), 
'(Uni/2), 1,1), ’(Z, 

V[1],V[2])»# 

The original transition equations five p(1,1,0) = 6p(0,0,0),.anticipating 
that p(1,3,0) may also occur, substitute p(l,3,0) = p(l,l,0) = 6p(0,0,0). 
Then use Qll to eliminate P(1,1,0). 

#'Q13^SUBSTITUTE(Q12, RIGHT(Qll), '(P(1,2,0)))# 

# ' Q14*~SIMPLIFY(SUBSTITUTE(Q13, '(THETA*P(0,0,0), THETA*P(0,0,0)) , 

*(P(1,3,0), P(1,1,0))))# 

Q14 contains only the unknown P(0,0,0) which can be determined from 
the condition G(l) = 1. 

#’Q15«-SIMPLIFY( SUBSTITUTE (RIGHT (Q14),l, f U) = 1)// 

The machine types out INDETERMINATE, indicating that 00 * 0 has been 
replaced by UNDEFINED. 

//EDISPLAY ( ’ Q15) # 

The operator LIMIT will be tried; this operator uses l’Hopital’s Rule. 
It is slow, and so it should not be used when substitution will suffice. 

#'Q16+-LIMIT(RIGHT(Q14) , 'U,l)// 

//EDISPLAY(' Q16) # 
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Q16 can be simplified by factoring and expansion. 

//'Q17-H?ACT0R0UT( EXPAND (Q16) , ’ (1((3*THETA + 2)*(1-2*THETA))))// 

//' Q17^LEFT(Q17) * EXPAND ( RIGHT ( Q17 ) ) // 

//EDISPLAY ( ' Q17 ) // 

Q17 is equal to 1 and can be solved for P(0,0,0). 

// ' Q18^S0LVE(Q17=l, 1 (P (0,0,0) ) ) // 

//EDISPLAY ( ' Q18) // 

The expression for P(0,0,0) shows that the SOLVE routine could be 
improved. This expression is now substituted into Q14 in order to produce 
the final expression for the generating function H. 

#' Q19*-SIMPLIFY(EVALUATE (Q14 , Q18) ) // 

Taking a census of Q19 shows that it is probably too large to display 
without being split. 

//CENSUS (Q18)// 

//EPRINT CLAST)// 1438 

In the thesis this expression was evaluated numerically to form a 
table of values. This probem would provide a basis for further work in 
automatic simplification. 




The preceding problem solutions show that the current system can 
be used for work on existing problems. No one part of the system is 
particularly weak, but there are many interesting possibilities for 
improvement. 



